The invention relates to image processing and in particular rendering of three-dimensional (3D) data, i.e. volume data. In particular, the invention relates to slab multi-planar reformatting rendering of volume data.
The process of calculating two-dimensional (2D) images of 3D objects is often referred to as volume rendering. Volume rendering finds applications in many fields. One such field is the rendering of medical volume data resulting, for example, from the scanning of the human or animal body using computed tomography (CT) and other X-ray scanners, nuclear magnetic resonance scanners and ultrasound scanners, to name but a few examples.
The volume data generated by modern scanning equipment can be very detailed and complex to interpret. A physician may wish to render the data using different view directions and from different positions with respect to the scanned object in order to be able to analyze the scanned object and to detect, for example, abnormalities.
Volume rendering techniques, such as applied to slab multi-planar reformatting (MPR) (sometimes referred to as MPR with thickness, or thick MPR), often lead to undesirable artifacts appearing in resulting 2D images. These artifacts are visually distracting and can hinder the interpretation of the images. In some situations the artifacts could be mistaken for real features of the volume data or in other cases could obscure real features of the data. Artifacts can also have a deleterious effect on subsequent image processing. For example, the accuracy of edge-detection algorithms is often very sensitive to the presence of such image artifacts.
Volume data comprises a plurality of voxels arranged in a 3D grid. Each voxel has a voxel value associated with it. The voxel values represent measurements of a physical parameter. For example, in the case of CT scans, the voxel values represent the opacity of those voxels to X-rays, i.e. their X-ray stopping power. X-ray stopping power is measured in Hounsfield units (HUs) which is closely correlated with density (mass per unit volume).
The voxels of volume data acquired by medical imaging devices are in most cases acquired on a Cartesian grid, i.e. the data points are aligned in three orthogonal axes I, J and K. In some special cases, the K axis may not be orthogonal to the I and J axes to take account of gantry tilt or slew. Moreover, the axes are conventionally ascribed a common origin at one corner of the volume. However, it will be appreciated that this choice of origin is arbitrary. These axes define a volume space. A volume-space coordinate system is used to identify the location of each voxel in the volume space. The volume-space coordinate system has unit (or basis) vectors i, j and k which are aligned with respective ones of the orthogonal axes I, J and K. The unit vectors i, j and k are defined such that the voxels are of unit length along each of the axes in volume space. That is to say, the separation between voxels (i.e. the distance between their centers) along each axis is unity.
The 2D image-data comprise a plurality of image pixels arranged in a 2D grid. Although the image itself is 2D, it is helpful to define a 3D view space containing the image. View space is defined by three orthogonal axes X, Y, Z having a common origin at one corner of the image. Again, the choice of origin is arbitrary. The X- and Y-axes are in the plane of the image (the image plane) and are aligned with the 2D grid of image pixels. The Z-axis is aligned parallel with the view axis (i.e. perpendicular to the image plane). A view-space coordinate system is defined to identify the location of each voxel and each image pixel in view space. The view-space coordinate system has unit, or basis, vectors x and y in the image plane and z along the view direction. The unit vectors x and y are defined such that the image pixels are of unit length along each of the axes in view space.
Images may be generated from the volume data using a conventional slab multi-planar reformatting (MPR) technique. In this technique, MPR data are generated by taking a coordinate in view space, transforming the coordinate into volume space, and resampling the volume data using some form of interpolation to generate a new MPR data value for the discrete view space coordinate. An MPR slice is formed by carrying this out for a plurality of {x, y} coordinates at a fixed z value. If this is repeated for multiple values of z, then multiple MPR slices are determined and these slices can be projected to form an MPR slab. The MPR slab thus comprises a series of MPR slices which are aligned parallel to the image plane and disposed at different positions along the Z-axis.
A 2D image for display is formed by projecting (collapsing) the MPR slab along the view direction onto the image plane. This is done according to a projection algorithm. The projection algorithm used in any particular case will depend on the desired appearance of the final image. Three different projection algorithms are commonly used. A first common projection algorithm is based on determining for each image pixel the maximum voxel value seen in the MPR slab along the Z-axis for the XY-coordinate corresponding to that image pixel. This is known as maximum intensity projection (MIP). Maximum intensity projection is a type of ray casting. In effect, for each pixel in the image, an imaginary ray is cast through the volume data parallel to the view direction. The image data for each pixel is then taken to be the maximum voxel value encountered by the ray as it traverses the MPR slab. Another common projection algorithm known as minimum intensity projection (MinIP) uses the minimum voxel value encountered by rays traversing the MPR slab for the image data instead of the maximum. A third type of projection mode used in slab MPR is Average Intensity Projection (AveIP) in which the voxel data values sampled from the portion of the ray traversing the slab are averaged to produce their collective value.
It will be appreciated that in some cases, only voxels having a voxel value in a selected range, or “window” will be of interest. For example, to reveal soft tissue on a CT scan, only voxel values in the range −200 to 500 HU may be of interest. To achieve such a view, a maximum or minimum intensity projection MPR image is typically calculated as described above, and subsequently the image is post-processed to enhance the contrast of voxel values in the desired range and suppress contrast outside that range.
Generally, the view axis will not coincide with one of the volume space axes, i.e. the image plane will have an arbitrary tilt angle with respect to any two of the volume space axes. Consequently sample points along a ray which is cast parallel to the view axis will not coincide with voxel coordinates, but rather each sample point along the ray will lie at arbitrary distances from a plurality of nearby voxel coordinates. Moreover, even if the view direction was along a volume space axis, in the general case each ray would not form a line which cut through a voxel coordinate in each slice of the volume data.
Because of this, when a view axis and a slab orthogonal to the view axis are defined, rays are generated which cut through the slab parallel to the view direction and spaced apart in the plane of the slab in a regular square grid. The rendering application then samples points along the portion of each ray that lies inside the slab, and for each sample point calculates a voxel value based on voxel values of the volume data that lie nearby. For example, for a ray passing almost equidistantly between adjacent voxel coordinates or centers, the average of the eight nearest voxel values, corresponding to the corners of a cube, might be used to arrive at the voxel value for the sample point. In general, a voxel value used for a particular ray passing through the MPR slice will be interpolated from the surrounding voxel values using a weighting based on the distance between the ray and the voxel centers. For example, a tri-linear interpolation between the eight surrounding voxels in the MPR slice is often used. Common interpolation methods include nearest neighbor (1 voxel), tri-linear interpolation (8 voxels), and tri-cubic interpolation (64 voxels).
A common and important type of slab MPR is so-called cine MPR in which the slab is gradually “moved” along the view axis in small increments, so that the user sees a succession of images or frames projected from the volume data, where each frame relates to a projection through a slab occupying a location incrementally different from the previous slab. In use, a radiologist will “cine” the slab back and forth interactively. In some cases, the succession of frames can be stored and then “played back” as a movie. For example, in a CT scan of a human head, a radiologist might take several seconds to traverse the head from top to bottom or back to front, with the traverse involving several hundred frames.
In cine MPR using MIP, MinIP or another non-linear projection function, the presence of noise in the volume data set will result in noise in the output image that changes from frame to frame along the view axis. As a practitioner cines a slab through the volume data set, the noise will cause a shimmering effect. This shimmering effect is distracting and also gives the impression of poor product quality. AveIP is less susceptible to this problem, because the noise is averaged out. Nevertheless AveIP is not entirely free from the shimmering problem, so the process may also be applied with AveIP.
As described above, a typical technique for generating such image frames is casting a ray through the volume data set for each pixel in an image and sampling the volume data sets at a number of discrete points along each respective ray. If using MIP, the maximum voxel value or interpolated value along the ray within the slab is then selected. The selected maximum sample or interpolated value within the slab is taken to be the pixel value associated with its respective ray.
FIG. 1 illustrates schematically a typical technique for producing a frame for cineing volume data. The volume data set 100 represents a volume space having a grid 102. The volume-space is illustrated in two-dimensions and comprises I and J axes, but also includes a third dimension along an axis K, which is perpendicular to the plane of the paper and therefore not visible. The voxel values for the volume data set are taken to be in the centre of each of the smaller squares and have a common origin at one corner of the volume data as illustrated by the dashed cross 104 in the lower left square of the volume space. Only 36 voxel values are illustrated, but it will be appreciated that many more values are typically used, but fewer have been used for illustrative purposes.
An image plane 106 is illustrated in FIG. 1 lying at an angle to the volume data set, i.e. tilted. In other words the image plane is not parallel to the volume data set. As described above, a number of rays 108 are cast from the image plane 106 into the volume data set 100 in a selected view direction parallel to a view axis Z. Each of the rays corresponds to a pixel at the image plane for forming an image. The image plane illustrated is 1-dimensional, for example along an axis X, but in practice the same process will be repeated in a second direction Y out of the plane of the paper. It is noted that in FIG. 1, and later figures of similar format, the K and Y axes are parallel. In reality, the K and Y axes have an arbitrary relationship and will not be aligned. However, to illustrate the principles clearly in a two-dimensional representation, the K and Y axes have been aligned.
A slab 114 is located in the volume-space 100 which is parallel to the image plane 106 and extends along the view axis Z. The slab is a 3D volume that also extends into the plane of the paper. The slab 114 is represented by a dashed rectangle. A number of discrete sample points 116 are selected that lie on each of the rays 108 and within the slab 114, which are illustrated as circles with a regular spacing that corresponds to the spacing of the rays 108. Before the sample points are collapsed or projected onto the image plane to form a 2D image, the values at each of the sample points is determined. For example, sample point 116 is illustrated as lying, by chance, near to the location of a voxel of the experimentally obtained data (i.e., near the centre of the square of the grids lines 102). The sample point 116 can therefore take for its voxel value the voxel value of the experimentally obtained voxel, or some value close to it. More generally, a sample point will not lie at or near to an experimentally obtained voxel location in which case an interpolation method will need to be applied to determine a suitable voxel value for the sample point from the voxel values of nearby experimentally obtained voxels. Sample point 118 depicts in two dimensions a more general case and lies approximately equidistant from four voxel locations. An interpolation method, such as tri-linear interpolation, might be used to obtain the value at the sample point 118. Once the values are obtained for each of the sample points, the sample point can be projected onto the image plane using MIP, for example, to obtain the values for each of the pixels.
FIG. 2 illustrates schematically an incremental movement of the slab from a first location 114 to a second location 120 as occurs during cineing. The incremental nature of the slab movement relative to the slab thickness is such that the majority of the “new” slab lies within the volume of the “old” slab, i.e. the two slabs 114 are 120 are largely overlapping. A new set of MPR planes is computed and projected for the slab 120 in the new position. The sample points of the old slab are shown as circles with a dashed line and the sample points of the new slab are shown circles with a solid line. For example, sample point 116 in the old slab corresponds to a shifted sample point 122 in the new slab. The voxel values in an experimentally obtained volume data set will of course include noise contributions. Because of noise, interpolated voxel values at sample points that are close together, such as sample points 116 and 122, can be significantly different. In the case that MIP is being used and the sample points 116 and 122 lie at the maximum voxel values for the ray in the slab, this will translate into significantly different brightnesses being assigned to the ray from one frame to the next. This is what gives rise to the shimmering effect. In other words, because the MIP result is dominated by a single sample point, even sub-pixel changes to sample point locations can produce significant intensity variations in the final intermediate image. The same is true for MinIP.
FIG. 3 shows a graph of values of sample points for the slabs shown in FIG. 2, along a single ray cast through the volume. The graph shows the sample values in arbitrary units against position along the ray, which is also in arbitrary units. The squares represent the sample points along the ray for the slab 114 at the first location and the triangles represent the sample points along the ray in the slab at the second location 120. Each of the square points and triangle points are spaced apart by 0.25 of a pixel from respective square and triangle sample points.
In FIG. 3, the solid vertical lines 134 and 136 show the boundaries of a slab of thickness τ in the first location and the dashed vertical lines 138 and 140 show the slab moved slightly to the right in a second location, where the size of the slab's movement is very small in relation to the thickness of the slab. However, despite the very small size of the slab movement, the magnitude of the maximum sample value 142 at the first slab position is significantly different from the magnitude of the maximum sample value 144 at the second slab location. As a result, when the user cines the slab through the volume, bright pixels in the MIP image appear to shimmer as the chosen sample points move closer or further from the theoretical maximum. A similar problem arises with MinIP, although in the illustrated example the magnitudes of the respective minimum intensity sample values, which are the two sample points contained in the ring 130, are very similar.